Microcanonical rates, gap times, 
and phase space dividing surfaces 

Gregory S. Ezraj^] 
Department of Chemistry and Chemical Biology 
Baker Laboratory 
Cornell University 
Ithaca, NY 14853 
USA 

Holger Waalken^ 

University of Croningen 
Institute of Mathematics and Computing Science 
P.O. Box 407 
9700 AK Croningen 
The Netherlands 

Stephen Wigging 

School of Mathematics 
University of Bristol 
Bristol BS8 ITW 
United Kingdom 
(Dated: January 20, 2009) 



1 



Abstract 

The general approach to classical unimolecular reaction rates due to Thiele is revisited in light 
of recent advances in the phase space formulation of transition state theory for multidimensional 
systems. Key concepts such as the phase space dividing surface separating reactants from products, 
the average gap time, and the volume of phase space associated with reactive trajectories, are both 
rigorously defined and readily computed within the phase space approach. We analyze in detail the 
gap time distribution and associated reactant lifetime distribution for the isomerization reaction 
HCN ^ CNH, previously studied using the methods of phase space transition state theory. Both 
algebraic (power law) and exponential decay regimes have been identified. Statistical estimates 
of the isomerization rate are compared with the numerically determined decay rate. Correcting 
the RRKM estimate to account for the measure of the reactant phase space region occupied by 
trapped trajectories results in a drastic overestimate of the isomerization rate. Compensating but 
as yet poorly understood trapping mechanisms in the reactant region serve to slow the escape rate 
sufficiently that the uncorrected RRKM estimate turns out to be reasonably accurate, at least at 
the particular energy studied. Examination of the decay properties of subsensembles of trajectories 
that exit the HCN well through either of 2 available symmetry related product channels shows that 
the complete trajectory ensemble effectively attains the full symmetry of the system phase space 
on a short timescale t < 0.5 ps, after which the product branching ratio is 1:1, the "statistical" 
value. At intermediate times, this statistical product ratio is accompanied by nonexponential 
(nonstatistical) decay. We point out close parallels between the dynamical behavior inferred from 
the gap time distribution for HCN and nonstatistical behavior recently identified in reactions of 
some organic molecules. 

PACS numbers: 05.45.-a, 82.20.-w, 82.20.Db, 82.30.Qt 
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I. INTRODUCTION 



The theory of unimolecular reaction rates, both for dissociative and isomerization pro- 
cesses, has been of great interest for nearly a century^. Following the original formulation 
of the statistical RRK and RRKM approaches to calculation of reaction rate^^'^'^'^, there 
has been a vast amount of activity, some of which has been described in several genera- 
tions of texbookj™^ i^ l i^ l i-^ l i^ l i-'^ l i« l ^^ l ^^ l ^™^'^^'^-^'^^ ] For concise overviews of the historical 
development of unimolecular rate theory, see ForstP^, Chapter 1, also ref. |25l 

To provide the context for the discussion in the present paper, we highlight some impor- 
tant and relevant contributions to the subject. (The selection of works cited is necessarily 
quite limited.) The development of Slater's "new" dynamical approach to classical uni- 
molecular dissociation rates^''^^ was followed by Thiele's general formulation of the problem, 
which emphasized the dynamical significance of the distribution of gap times (see below 
As we show below, Thiele's work has proved to be remarkably prescient in terms of its 
identification of the appropriate phase space structures involved in unimolecular reaction 
dynamics. Following Thiele's work, a large number of computational (trajectory) investi- 
gations of lifetime/gap time distributions have been undertaken (see, for example, results 
discussed iiJEMEll^ Dumont and Brumei'^^'^ and DumonlplEZEsEl j^g^yg analysed unimolec- 
ular reaction rates in terms of the gap time distribution, while related work has been done 
on the so-called classical spectral theoremp5ESI_ 

A fundamental assumption of the statistical theory of (classical) unimolecular decay is 
that intramolecular vibrational energy distribution (IVR) occurs on a timescale much faster 
than that for reactionf^^'^. Renewed interest in the problem therefore naturally stemmed 
from investigations of the properties of molecular vibrational dynamics, modelled as non- 
linearly coupled anharmonic oscillators, in light of the KAM theorem and the apparent 
existence of a threshold energy for onset of global chaoj^^ ^ * ^^ * ^° l The central role of deter- 
ministic chaos itself in determining the validity of statistical approaches to reaction rates, as 
well as possible quantum manifestations of classical nonintegrable behavioi"^, has received 
much attentiorP^, most recently in the context of quantum controP^HS ^j^g possibility 
of mode-specific chemistry also stimulated much work on the relation between intramolec- 
ular dynamics and non-statistical reaction dynamics^^. Several examples of non-RRKM 
effects in thermal reactions of medium-sized organic molecules have been found in recent 
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yeard^SESIiaiH!^ while a detailed discussion of possible non-RRKM effects in one particular 
reaction, the unimolecular dissociation of H^, has been givenH^l. 

Advances in the classical theory of chemical reaction rates have been closely linked 
to improvements in our understanding of the phase space structure of Hamiltonian 
systemd2^'^^'S2I13_ Both conceptually and in practice it is important to distinguish between 
systems with 2 degrees of freedom (DoF) and multimode (A^ > 3 DoF) systems. 

For 2 DoF systems described by standard kinetic plus potential Hamiltonians, periodic 
orbit dividing surfaces (PODS; that is, dividing surfaces in configuration space obtained by 
projection of a periodic orbit from phase space) were shown to have the property of minimal 
flux (as a consequence of the principle of stationary action) and to be (locally) surfaces of 
no return (the velocity vector is nowhere tangent to the configuration space projection of 
the periodic orbit )^^'^^'^. In this case the PODS therefore provide a rigorous realization of 
the concept of the transition statd^^'^. (PODS have also been defined for 2 DoF systems 
lacking time-reversal symmetry^.) Correlations between features of the classical phase 
space structure and the behavior of the computed reactive fiuxP^'^ were explored relatively 
early on for 2 DoF systems by DeLeon and Bernd^^ESI^ Important theoretical advances were 
subsequently made relating non-statistical behavior to molecular phase space structur^^, 
in particular the existence of intramolecular bottlenecks to energy transfej ^^ * ^^ * ^^ !, broken 
separatriced^S'SS^, and reactive islands and cylinder j^^ ' ^^ ' ^^ ' ^^^ ' ^^ ' ^^ ' ^^ l 

The multimode case remains less thoroughly explored and understood. Methods for 
defining approximate intramolecular bottlenecks and reactive dividing surfaces have been 
dg^igg jMMSQlHIIHaMIHll. Although the Arnold web of resonances provides a useful frame- 
work for mapping and analyzing evolution of phase points in near-integrable multimode 
systemd^l'^, so-called "Arnold diffusion'''^^'^ has become a convenient but often ill-defined 
catchall term employed to describe a variety of possibly distinct phase space transport 
mechanisms in A^ > 3 DoF system^. The possible role of "Arnold diffusion" in multimode 
molecular systems has been studied in iVRpSElIHHISQllil ^^^^ isomerizatioEp^'^. Local ran- 
dom matrix theories have served as the foundation for quantum theoretical treatments of 
isomerization reactions in large molecule d^^*^^ * ^^*^ -^*^, while attention has also been given to 
related scaling approache d^^ * ^°° * ^ -^ and to fractional kinetic4 ^°^ * ^°^ * ^°^ l 

An important advance was the realization that normally hyperbolic invariant manifolds 
(NHIMj^^SnnSl^ provide a natural and theoretically well-founded generalization of PODS to 
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the > 3 DoF case (cf. ref. ll07p . The stable and unstable manifolds associated with NHIMs 
define co dimension-one dividing surface^^ on the constant energy 

manifolcPSnei^ and so 

are possible candidates for reactive separatrices, while the NHIMs themselves define phase 
space transition state^^^^'^. Some early attempts were made to compute and visualize such 
manifolds in a 3 DoF system describing surface diffusion of atom^^ and for a 4D symplectic 
mapping modelling the dissociation of a van der Waals complex^^^EE] 

As discussed in more detail below, based on the notion of the NHIM and on the de- 
velopment of efficient algorithms for computing normal forms at saddles, there has been 
significant recent progress in the development and implementation of phase space transition 
state thcor }^^"^'^"'^'^^"'^^^'^^^'^^^'^^^'^^^'^^^'^^'^' (see ab j^^^'^^^'^^^'^^^'^^^'^^^'^^^'^^^""^ . 

In the present paper we consider the problem of defining and evaluating theoretical 
unimolecular reaction rates in light of the penetrating analyses of Thield^Zl^ Dumont and 
Brumei'^, and DeLeon and Bern#^'^, the classical spectral theorempSESl, and the re- 
cent developments in phase space transition state theory^ ^ l ^^^ * ^^"^ * ^^^ * ^^^ * ^^^ * ^^^ * ^^^ * ^^*^ ! men- 
tioned above. The particular reaction chosen for study is the isomerization HCN 
Q^;Hl ii(i | ii7 | i29 | i30 | i3i | i32 | Kj3 | i34 | i35 | i36 | i37 | ^j^j^j^ scvcral relevant theoretical quantities have 

recently been computed for a (classical) model of the HCN molecule at fixed energji'^I^'^. 

We are concerned with the rate of reaction, at fixed energy, for a system described by 
a time-independent, n degree-of-freedom (DOF) classical Hamiltonian. One measure of the 
rate at which trajectories leave a region of the energy surface is given by the (magnitude of 
the) fiux of trajectories leaving that region (units of energy surface volume/time) divided by 
the energy surface volume of initial conditions in that region corresponding to trajectories 
that will eventually leave the region. This rate is just the inverse of the mean passage or gap 
tim^2Zl29!_ Hence, to compute this rate at fixed energy one must first (1) define the region 
of reactants, then (2) compute the fiux of trajectories exiting this region and, finally, (3) 
compute the volume of the energy surface corresponding to initial conditions of trajectories 
that leave the reactant region. 

The importance of the gap time distribution was emphasized by Thield^, who ex- 
plicitly invoked the concept of a phase space dividing surface separating reactants and 
product j^^^ * ^^^ l While the use of dividing surfaces (transition states) defined in configuration 
space is quite familiar in the field of reaction rate theorj^SEIIlQ]^ carrying out these steps in 
phase space, as opposed to configuration space, remains less familiar in practicd^^'^. In the 
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present paper we show how the expression for the microcanonical rate of reaction described 
above can be evaluated using the phase space approach to reaction dynamics developed 
in a recent series of paper d^"^ ' ^^" ' ^^^ ' "^ ' "^ ' ^^^ ' ^^^ ' "^ ' ^'" '. Moreover, we analyze in detail the 
properties of the gap time distribution previously obtained for HCN isomerization using the 
phase space reaction rate theorjW^. 

The stucture of the paper is as follows: in Section [IT] we outline an approach to 
the classical theory of unimolecular reaction rates based upon the general formulation 
of Thiel^^. We demonstrate that various dynamical quantities related to the unimolec- 
ular reaction rate defined within Thiele's approach can be rigorously defined and com- 
puted within the recently developed theoretical framework for phase space transition state 
l;]^QQ];.-ylio 9 | iio | ii4 | ii5 | ii6 | ii7 | ii8 | ii9 | i20 | gpgcial emphasis is placed on the properties of the gap 

time distribution, and the relation between the inverse gap time and the rate of reaction. In 



Section III we consider the dynamics of the HCN isomerization reaction. We analyze pre- 
viously computed classical trajectory result^nSEZl jj^ light of the discussion given in Section 
[IT| Both gap time and reactant lifetime distributions are analyzed, and distinct temporal 
regimes identified, with power law (algebraic) decay seen at intermediate times and expo- 
nential decay at long times. Moreover, we find a statistical (1:1) product ratio for both 
algebraic and exponetial decay regimes. Comparisons are made between computed values 
of the decay rate and statistical estimates of the isomerization rate. Our results show that 
analysis of gap time distributions computed using phase space dividing surfaces having min- 
imal flux and no-recrossing properties yields a great deal of information about the system 



dynamics. Section IV concludes. 
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II. GENERAL APPROACH TO UNIMOLECULAR REACTION RATES 



In this section we present a general formulation of unimolecular reaction rates based 
upon that originally given by Thield^. In their general form the rate expressions derived by 
Thiele explicitly invoke the existence of a phase space dividing surface separating reactants 
and products; such surfaces, although discussed by Wignei'^ at the inception of the theory 
of chemical reaction rates (see alscP^SHini^ ^ have only recently become amenable to direct 
computation via the use of normal form approache^^ ^" ' ^^^ ' "^ ' ^^^ ' ^^^ ' ^^^ ' "^ ' ^"" ' ^^^ ' ^^^ ' ^ ^ 

A. Phase space dividing surfaces: definition and properties 

We are concerned with the rate of unimolecular reactions, either dissociation or isomer- 
ization, at fixed energy, for a system described by a time-independent, n degree-of-freedom 
(DOF) classical Hamiltonian. Points in the 2n-dimensional system phase space Ai = M^" 
are denoted x & Ai. The system Hamiltonian is H{x), and the {2n — 1) dimensional energy 
shell at energy E, H{x) = E, is denoted E^; C Ai. The corresponding microcanonical phase 
space density is 6{E — H{x)), and the associated density of states for the complete energy 
shell at energy E is 

p{E)= [ dx 6{E - H{x)). (2.1) 
Jm 

The first step in the analysis is to define the region of the energy surface corresponding to 
the reactant of interest. For complex systems there may be many such regions of interest!^, 
but for the moment we focus on only one such region, and the rate at which trajectories leave 
that particular region. We treat either unimolecular dissociations or isomerizations in which 
molecules are removed from consideration immediately after passing from the reactant region 
of phase space, and are not therefore directly concerned here with reactive flux correlation 
functions or associated relaxation kinetic ^^ * ^^ * ^" * ^^ * ^ -^. 

A standard approach to defining the reactant region is to consider the potential energy 
function-a configuration space based approach?^. In the simplest case the reactant region 
is associated with a single local minimum, and trajectories exit/enter the reactant region by 
passing over energetically accessible saddle points of the potential. Although there could be 
several saddles controlling rates of exit from the potential well, as mentioned we will only 
consider the case of a single saddle. We comment on the case of multiple saddles below. 
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However, dynamics occurs in phase space (i.e., on the energy surface S^; C for the 
situation we are considering) while the picture described above is concerned with the pro- 
jection of trajectories from phase space to configuration space. For multidimensional sys- 
tems such as polyatomic molecules (A^ > 3 DoF), it is in general no longer possible to 
define or compute a dividing surface with desirable dynamical attributes such as the no- 
recrossing property by working in configuration space alone, and a phase space perspective 
is ncccssar} ^^"^'^^"'^^^'^^^'^^^'^^'''^^^'^^^'^^"'^^^'^^^'^^^' 

While the phase space approach to reaction dynamics reviewed briefly below does not 
give a complete solution to the problem (questions of transition state bifurcations and non- 
integrability remain to be addressed, for example), it does provide an important component, 
which, when coupled with knowledge of certain properties of the system Hamiltonian func- 
tion, does provide the necessary information to identify the reactant region in some situations 
of interest. Saddle points of the potential energy surface still play a role - in particular, 
rank one saddles. For Hamiltonian functions that are the sum of the kinetic energy and the 
potential energy a rank one saddle point of the potential energy function is manifested as 
an equilibrium point of saddle type for the 2n dimensional Hamilton's equation^^. More 
precisely, it is an equilibrium point of saddle-center-. . .-center stability type, meaning that 
the 2n x 2n matrix associated with the linearization of Hamilton's equations about this 
equilibrium point has 2n eigenvalues of the form ±A, izcufc, = 2, . . . , n (A, tUfc > 0). The 
significance of saddle points of this type for Hamilton's equations is that, for a range of 
energies above that of the saddle (we explicitly discuss the range of energies later on), the 
energy surfaces have the bottleneck property in a phase space neighborhood near the saddle, 
i.e. the 2?2 — 1 dimensional energy surface locally has the geometrical structure of the prod- 
uct of a 2n — 2 dimensional sphere and an interval, 5*^""^ x I. The bottleneck property is 
significant because in the vicinity of the bottleneck we are able to construct a dividing sur- 
face, DS(i?), (where 'E' denotes the energy dependence) with very desirable properties. For 
each energy in this range above the saddle DS(i?) locally "disconnects" the energy surface 
into two, disjoint pieces with the consequence that the only way to pass from one piece of 
the energy surface to the other is to cross DS(-E'). The dividing surface has the geometrical 
structure of a 2n — 2 dimensional sphere, S"^""^, which is divided into two 2n — 2 dimensional 
hemispheres, denoted DSin(i?) and DSout(-£') that are joined at an equator, which is a 2ra — 3 
dimensional sphere, S'^"^'^. The hemisphere DSin(-E') corresponds to initial conditions of 
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trajectories that enter the reaction region while DSoutl-E") corresponds to initial conditions 
of trajectories that exit the reaction region, both by passing through the bottleneck in the 
energy surface. The equator S*^""^ is an invariant manifold of saddle stability type, a so- 
called normally hyperbolic invariant manifold (NHIM)^. The NHIM acts as the "anchor" 
for this entire construction and is of great physical significance: it is the actual "saddle" in 
phase space identified as the "activated complex" of reaction rate dynamic . Our 
focus here is on microcanonical rates, and it has been shown that DSin(-E) and DSout(-E') 
have the essential no-recrossing property and that the flux across them is minimal^. We 
denote the directional flux across these hemispheres by (t>in{E) and (f)o^t{E), respectively, and 
note that (f)m{E) + 0out(-E) = 0. For our purposes we only need the magnitude of the flux, 
and so set |0in(-E)| = \(f)out{E)\ = (p{E). Most significantly, the hemisphere DSin(-E') is the 
correct surface across which to compute the "exact" flux into the reaction region. 

B. Phase space volumes and gap times 

The disjoint regions of phase space corresponding to species A (reactant) and B (product) 
will be denoted A^a and A^b, respectively'Hl, We assume that all coordinates and momenta 
are bounded on the reactant energy shell S^; fl M.a, and that it is possible to define a 
boundary (dividing surface) DS in phase space separating species A and B. As discussed 
above, the DS can be rigorously defined to be locally a surface of no return (transition state). 
The microcanonical density of states for reactant species A is 

pj,{E)= [ dx5{E-H{x)) (2.2) 

J Ma 

with a corresponding expression for the density of states Pb{E) for product B for the case 
of compact product energy shell tV^b- 

Provided that the flow is everywhere transverse to DSin, out(-E'), those phase points in the 
reactant region A^a that lie on crossing trajectorie^^^'^ (i.e., that will react, and so are 
"interesting" in Slater's terminology^) can be specified uniquely by coordinates {q,p,ilj), 
where {q,p) G DSin(-E') is a point on DSin(-E'), the incoming half of the DS, specified by 
2{n — 1) coordinates {q,p), and ip is a time variable. (Dividing surfaces constructed by 
the normal form algorithm are guaranteed to be transverse to the vector field, except at 
the NHIM, where the vector field is tangent! ^°^ * ^^° l) The point x{q,p,i/j) is reached by 
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propagating the initial condition {q,p) E DSin(-E') forward for time ip (see Figure [T]). As 
all initial conditions on DSin(-E) (apart from a set of trajectories of measure zero lying on 
stable manifolds) will leave the reactant region in finite time by crossing DSout (-£'), for each 
(^)P) ^ DSin(-E) we can define the gap time s = s{q,p), which is the time it takes for the 
incoming trajectory to traverse the reactant region. That is, x{q, p,ip = s{q, p)) G DSout(-E). 
For the phase point x{q,p,ilj), we therefore have < ■?/' < s{q,p). 

The coordinate transformation x — > {E, tp, q,p) is canonica P^ * ^^^ * ^"^^ -*^, so that the phase 
space volume element is 

d^''x = dEdtljda (2.3) 

with da = d"~^gd"~^p an element of 2n — 2 dimensional area on the DS. 

As defined above, the magnitude (p{E) of the flux through dividing surface DS(£^) at 
energy E is given by 



0(E) 



(2.4) 



da 

'DSi„(£;) 

where the element of area da is precisely the restriction to DS of the appropriate flux 
(2n — 2)-form ct;*^"~^Y(n — 1)! corresponding to the Hamiltonian vector field associated with 
^^^•j[IIIIIl5lT46lii3_ rpj^g reactant phase space volume occupied by points initiated on the 
dividing surface DSin with energies between E and E + dE is therefor d^^ * ^^ * ^^ * ^^^ * ^^^ * ^^^ * ^ ^ 



dE f da f di) =dE ! da s (2.5a) 

= dE(l){E)s (2.5b) 

where the mean gap time s is defined as 

1 r 

da s (2.6) 



1 



0(E) 



DSinCS) 



and is a function of energy E. The reactant density of states p^{E) associated with crossing 
trajectories only (those trajectories that enter and exit the reactant region|S2; see below) is 
then 

pI{E)=cP{E)-s (2.7) 



where the superscript C indicates the restriction to crossing trajectories. The result (2.7) is 
essentially the content of the so-called classical spectral theoreni ^^ * ^^ * ^^^ * ^^^ * ^^^ * ^^^ l 

If all points in the reactant region of phase space eventually react (that is, all points lie 
on crossing trajectorie d^^ * ^^ *) then Px{E) = pa{E), the full reactant phase space density of 

10 



states. Apart from a set of measure zero, all phase points x G M.a can be classified as either 
trapped (T) or crossing (C)^^. (Further discussion of this division of the reactant phase 
space in terms of the Poincare recurrence theorem is given in Appendix |Aj) A phase point 
in the trapped region M.\ never crosses the DS, so that the associated trajectory does not 
contribute to the reactive flux. Phase points in the crossing region do however eventually 
cross the dividing surface, and so lie on trajectories that contribute to the reactive flux. In 
general, however, as a consequence of the existence of trapped trajectories (either trajectories 
on invariant trapped n-tori^^'^ or trajectories asymptotic to other invariant objects of zero 
measure), we have the inequality'^^'^ 

pI{E)<pp.{E). (2.8) 

If p^{E) < pa{E), then it is in principle necessary to introduce corrections to statistical 
estimates of reaction rate j^^ * ^^ * ^" * ^^^! Numerical results for p'^{E) and p{E) for the HCN 
molecule are discussed belo'wW^*^. 

C. Gap time and reactant lifetime distributions 

Of central interest is the gap time distribution, V{s;E): the probability that a phase 
point on DSin(-E') at energy E has a gap time between s and s + ds is equal to V{s; E)ds. 
An important idealized gap distribution is the random, exponential distribution 

V{s;E) =k{E)e-^^'^^' (2.9) 

characterized by a single decay constant k (where k depends on energy E), with correspond- 
ing mean gap time s = k~^. 

The lifetime (time to cross the dividing surface DSout(-E')) of phase point x{q,p,ip) is 
t = s(q,p) — ip (cf. Fig. [T|d). The volume of reactant phase space occupied by trajectories 
having lifetimes t > t' at energy E is then 

r+oo 

Vol{t>t';E)=(f){E) / ds {s -t')V{s;E) (2.10) 

Jt' 

so that the corresponding probability of an interesting phase point in the reactant region 
having a lifetime t > t' is obtained by dividing this volume by the total volume occupied by 
points on crossing trajectories, (p{E)s, 

PToh(t>t';E) = - / dsV(s;E)(s -t'). (2.11) 

s Jt' 
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The corresponding reactant lifetime distribution function P(t; E) at energy E is therefore 

P(t; ^) = - A Prob(t > t'; E) ^ (2.12a) 
= - / dsVis-E) (2.12b) 

s Jt 

where the fraction of interesting (reactive) phase points having hfetimes between t and t + dt 



is P(t; E)dt. It is straightforward to verify that the hfetime distribution (2.12) is normahzed. 



Note that an exponential gap distribution (|2.9|) imphes that the reactant hfetime distribu- 
tion 



E) is also exponential; both gap and lifetime distributions for realistic molecular 
potentials have been of great interest since the earliest days of trajectory simulations of 
unimolecular decaj^^^*^. 



We emphasize that the rigorous relation (2.12) between the gap time distribution and 



the reactant lifetime distribution follows quite straightforwardly from our continuous time 
formulation of the problem using the canonical transformation of phase space variables eq. 



(2.3) and the properties of the dividing surfaces DS (cf. ref. [29]). 

A concise derivation of the delay differential equation for the Delayed Lifetime Gap 
ModeP^ is presented in Appendix [b] 



D. Reaction rates 



We start with the (classical) expression for the rate k{T) of a collisionally activated uni- 
molecular process at temperature T derived by ThieldS'^SEZl^ Using the notation established 
above, the rate coefficient k{T) is 

-OO P 

/ da [1- e~^'] (2.13a) 



k{T) = ^ I dEe- 

JEo ^DSin(S) 



-I3E 



Here, Z\{T) is the reactant partition function 



1 P+CO 

- dE e-^^ (f){E) [1 - ^] 

'A J En 



(2.13b) 



Za = JdE e-^V(^), 



(2.14) 



(5 = l/k-sT, uj is the effective collision rate per molecule, Eq is the threshold energy for 



reaction, and the overline in eq. (2.13b) denotes an average over the dividing surface DSin(-E). 



The physical interpretation of expression (2.13 ) is that the thermal reaction rate k{T) is given 
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by the average of the equihbrium activation rate times the probabihty that an activated 
phase point will react (that is, cross the dividing surface DSout) before it suffers a (strong) 
deactivating collision. 



We note that the rate expression (2.13) makes sense even though for larger energies E 
the dividing surface DSin(-E') with desired dynamical properties might no longer exist, as 
the contribution from high energies is damped away by the exponential Boltzmann factor. 
(Thiel^^ assumed the existence of a suitable dividing surface for all energies E.) 

The limiting expressions for k{T) obtained at high and low pressures are of interest. At 
high pressures (a; ^ oo) we have 



fcoo(T) = lim k{T) = — dE e-f^^ / da (2.15a) 

'^^^ Jeo JT)Sn,{E) 

1 r+oo 

= — / dEe-^^0(E) (2.15b) 

^A Jeo 
1 r+oo 

= — / dEpA(^)e-^^e™(^) (2-15c) 
Za Jeo 



where the quantity 



ifc™(i?) ^ (2.16) 

is the statistical (RRKM) microcanonical rate for the forward reaction (A B) at energy 
E, the ratio of the magnitude of the flux (f){E) through DSin(£') to the total reactant density 
of state^^^^. The rate koo{T) is then the canonical average of the microcanonical statistical 
rate kJ^^^{E). (The collision rate ut should not be so large that trajectories of systems 
crossing the dividing surface are significantly perturbed by collisions.) 
Clearly, if pa{E) = p2(^)> then 

1 

s 

the inverse mean gap time. In general, the inverse of the mean gap time is 



fc™(E) = ^ (2.17) 



^-^^^^ kf^^M ^2.18a) 



Pa 



RRKM 



Pa{E) 

p1{e) 



(2.18b) 



> kf^"™. (2.18c) 

The rate k^^^ can be interpreted as the statistical unimolecular reaction rate corrected 
for the volume of trapped trajectories in the reactant phase spac# ^^ * ^^ * ^° l The modified 
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statistical rate is therefore predicted to be greater than the RRKM rate, a prediction usually 
at odds with numerical simulations (cf. results presented below, also refs H7P^63ll7Up . 
The low pressure {u 0) limit of the rate is 

ko{T) = lim k{T) = — dE e'^^ / da s (2.19a) 

= — / dEe-^^(P{E)s (2.19b) 

Jeo 
1 r+oo 

= ^ / dE pA(^)e-^^ A;™(E) us (2.19c) 
Za Jeo 

showing that the effective microcanonical rate is smaller than the RRKM statistical weight 
by a factor u;s ^ 1. That is, in the low pressure limit the reaction rate is proportional to 
the rate of collisional activation u, with each molecule taking a time s on average to react. 

Different choices of transition state location will in general result in different mean gap 
time^. We emphasize that exact and unambiguous calculation of the mean gap time 



(2.18) is possible given knowledge of the phase space geometrical structures that enable 
us to construct the reaction region, a dividing surface with minimal flux (hence the exact 
flux can be computed without integrating trajectories), and the reactive volume, i.e., the 
volume of the energy surface corresponding to interesting initial conditions x G A^a- The 
fundamental geometrical structures required to compute these quantities are the 2n — 2 
dimensional hemipheres DSin(-E') and DSout(-E') that control the entrance to the reaction 
region and exit from the reaction region, respectively. These geometrical structures are 
what we use to compute flux and we sample initial conditions on DSin(-E') and integrate 
them until they reach DSout(-E) in order to compute SDSi„(£;). ^ detailed algorithm has 
previously been given for computing the dividing surfaces and the flux across 
In these references, numerical tests were also described for determining the range of energies 
above the saddle for which the dividing surface "locally disconnects" the energy surface in 
the way described above (in particular, for energies sufficiently larger than the saddle the 
dynamics may not "feel" the influence of the saddle point at all and the energy surface could 
deform so much that it would make no sense to speak of disjoint regions of "reactants" and 
"products"'H3i 
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E. Multiple saddles 



So far, we have only considered a reaction region where access in and out of the region is 
controlled by a single saddle. We can similarly consider a reaction region (in phase space) 
where access in and out of the region is controlled by d saddles (actually, d saddle-center- . . . 
- center type equilibria). Associated with each saddle we can compute a2n — 2 dimensional 
dividing surface, DSi(E), that is divided into two 2n — 2 dimensional hemispheres DSi,in(-E) 
and DSi^out(-E') which have the same interpretation as above. The magnitude of the flux out 
of the reactive region is denoted by XliLi 4>i{E) and it is shown vt^^^^ihaX the volume of 
initial conditions in the energy surface that correspond to trajectories that leave the reaction 
region is given by XliLi ^DSi !„(£;) 0i(-£') 7 where sdsi the average time for trajectories 

starting on DSi4n(-E') to cross DSj,out(-E'), for any 1 < j < d. In this case the corrected total 
statistical escape rate is given by: 

fc™(E; d) = .^-^'^^^^l^^ ■ (2.20) 

An important case is that where all of the saddles are symmetric in the sense that 
HE) = HE), (E) = SDSj,i„(£;) 0j,in(-E), for all 1 < i,j < d then the corrected 



statistical rate (2.20) reduces to expression (2.18). 

The symmetric situation applies to the case of HCN isomerization considered irfUMHI 



Numerical results for HCN are further discussed in Sec. Ill of the present paper. 
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III. HCN ISOMERIZATION DYNAMICS 



The isomerization dynamics of HCN, HCN ^ CNH, has been widely studied, using both 
classical and quantum mechanics: see, for example, refs I116II117|129|13U|I131II132|133|1341 
11351)136111371 and references therein. In the calculations reported in refs lll6llll7l the potential 
energy surface of Murrell, Carter, and Halonen (MCH)^^° was used. For the MCH potential 
energy surface the saddle point is at energy -12.08 eV and the trajectory calculations in refs 
111611171 were all carried out at energy 0.2 eV above the saddle. 

The HCN molecule is modelled as a planar system with zero angular momentum, so 
that there are three DoF ^^^ * ^^ . In planar HCN there are two saddles, related by reflection 
symmetry, separating reactant (HCN) from product (CNH), with bond angle 7 = ±7* ~ 
±67° (see Figure [2]). The mean gap time is found to be s = 0.174 ps which corresponds to 
an isomerization rate of 0.14 x 10~^ a.uJ ^^^ * ^^ ^ (see below). A discussion of numerical aspects 
and efficiency of the calculations was also given in refs I116|I117[ 

A. Gap time and reactant lifetime distributions 

The phase space structures of interest at fixed energy E, namely the NHIMs and the di- 
viding surfaces separating reactant from product, are computed via an algorithmic procedure 
based on Poincare-Birkhoff normalization that is described in refs I110I120I The Poincare- 
Birkhoff normalization provides a nonlinear, symplectic transformation from the original, 
physical coordinates {q,p) to a new set of coordinates, the normal form coordinates {q,p), 
in terms of which the dynamics is "simple" . Morover, in these normal form coordinates the 
phase space structures governing reaction dynamics can be expressed in terms of explicit 
formulae, as described in refs lll0|ll20[ Their infiuence on the dynamics, in the normal form 
coordinates, is very easy to understand, and the geometrical structures can then be mapped 
back into the original coordinate system. 

The Poincare-Birkhoff normal form theory provides an algorithm to compute the sym- 
plectic tranformation T from physical coordinate to normal coordinates, 

T{q,p) = {q,p). (3.1) 

In a local neighbourhood C of the equilibrium point of interest, this trandformation "un- 
folds" the dynamics into a "reaction coordinate" and "bath modes" : expressing the system 
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Hamiltonian H in the new coordinates, {q,p), via 

HnF{q,P)=H{T-\q,p)), 



(3.2) 



gives i^NF in a simplified form. This "unfolding" into a reaction coordinate and bath modes 
is one way of understanding how we are able to construct the phase space structures, in the 
normal form coordinates, that govern the dynamics of reaction. The explicit expressions for 
the coordinate transformations, T{q,p) = {q,p) and T~^{q,p) = {q,p), between the normal 
form (NF) coordinates and the original coordinates provided by the normalization procedure 
are also essential, as they allow us to transform the phase space structures constructed in 
normal form coordinates back into the original physical coordinates. 

We consider an ensemble of trajectories with initial conditions sampled uniformly on one 



of the two dividing surfaces, DSi^in(-E) say, according to the measure da (cf. eq. (2.4)), so 
that all trajectories initially enter the reactant region of phase space via channel 1. Initial 
conditions on the dividing surface are obtained by uniformly sampling the 2n — 2 dimensional 
DS in normal form coordinates (g,p) with reaction coordinate variables (^i,pi) determined 
by energy conservation {E = —11.88 eV) together with the constraint qi = pi. The inverse 
of the normal form coordinate transformation is used to compute physical coordinates for 
initial phase points, and trajectories are propagated forward in time in physical coordinates. 
As a trajectory approaches either dividing surface DSj,out(-E'), a transformation is made to 
normal form coordinates, which allows accurate determination of the time at which the 
trajectory crosses the out DS. This is the gap time. 

The ensemble consists of 815871 trajectories; of these, 83599 ultimately exit through 
channel (DS) 1, while 732272 exit through channel 2. The branching ratio for the whole 
ensemble is then 8.76, very different from the statistical value unity dictated by symmetry. 
The mean gap time for the complete ensemble is 0.174 p4^. For trajectories reacting 
via channel 1, s = 0.714 ps, while s = 0.112 ps for those reacting via channel 2. Those 
trajectories that exit via the same transition state through which they entered therefore 
have a significantly larger average gap time. This makes physical sense, as such trajectories 
must have at least one turning point in the bending motion. 

The gap time distribution V{s; E) for the complete ensemble at constant E = —11.88 eV 
is shown in Figure |3^ (cf. Figure 4b of ref. I117p . Gap time distributions for subensembles 1 
and 2 are shown in Figures [3]d and [3}:;, respectively. In addition to the gap time distribution 
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itself, we also consider the cumulative distribution F{t), which is defined as the fraction of 
trajectories on the DS with gap times s > t, and is simply the product of the normalized 



reactant lifetime distribution function P(t; E) and the mean gap time s (cf. eq. (2.12)): 

/ + 00 
dsV{s;E) (3.3a) 

= sF{t;E). (3.3b) 



For the random gap time distribution (2.9), the cumulative gap time distribution is expo- 
nential, F{t) = e~^\ 

Cumulative gap time distributions for < t < 0.4 ps are shown for the whole ensemble 
in Fig. |4^, and for the two subensembles in Figs |4)d and |4]3, respectively. Reactant lifetime 
distributions for the whole ensemble over longer time intervals (0 < t < 25 ps) are presented 
in Figure |5] both log[F(t)] vs t and log[F(t)] vs logt plots are shown. Corresponding plots 
for the two subensembles are shown in Figs |6] and [7j 

We now discuss properties of the gap time and lifetime distributions on various physically 
relevant timescales: 

1. Very short times (t <^ s) 

By construction, the phase space DS used to separate reactants and products elim- 
inates local (short-time) recrossing ij^"^ ' ^^" ' ^^^ ' ^^^ ' ^^^ ' ^^^ ' ^^^ ' ^^^ ' ^ . Hence, there are no 
very short gaps. 

2. Short times (t ~ s). 

The gap time distribution for the complete ensemble shows "pulses" of reacting tra- 
jectories. Each pulse is associated with a bundle of trajectories that execute a certain 
number of oscillations in the reactant well before crossing one or the other Dg ii6 | ii7 j_ 
The first pulse is associated with trajectories exiting via channel 2, the second pulse 
with trajectories exiting via channel 1, and so on. Similar structure has been seen in 
the lifetime distribution computed for escape of Rydberg electrons in crossed field^^. 

The cumulative gap time distributions up to times ~ s (Fig. |4]) exhibit a structured 
and faster-than-exponential decay and so cannot readily be fitted to an exponential 
curve in order to obtain an effective decay rate for trajectories leaving the HCN well. 
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3. Intermediate times 10s '^t>s. 

On intermediate timescales, the peaks associated in the gap time distribution associ- 
ated with individual pulses begin to overlap. The reactant lifetime distribution exhibits 
algebraic (power law) decay, F(t) ~ 1/t", with a ~ 0.705 (Figure |5^). Such power 
law decay has been seen in other models for isomerizatiorP^^*^, and is in general 
associated with fractional kinetic^^. 

Attempts have been made to rationalize the existence of power law lifetime distribu- 
tions in such systems in terms of a hierarchical set of bottlenecks, perhaps associated 
with the Arnold wetP^^ presumed to exist in the vicinity of the minimum of the po- 
tential well. At this stage, however, a more quantitative explanation of the dynamical 
origins of algebraic decays such as those seen here at intermediate times remains an 
open problem. 

On the same timescale, both subensembles exhibit algebraic decay with essentially 
identical exponents a: 0.708 and 0.701, respectively (Figs [6^, [7^). Of course, if one 
starts with an ensemble of reactant phase points whose distribution possesses the 
symmetry of the phase space induced by the reflection symmetry of the potential, 
then equality of the exponents for the algebraic portion of the decay is expected. Our 
ensemble of initial conditions is nevertheless highly asymmetric. We comment on the 
observed exponent equality further below in our discussion of branching ratios. 

4. Long times t s. 

At longer times the lifetime distribution exhibits exponential decay, F{t) ~ e"*^*, 
with exponent k ~ 0.092 ps~^ (the decay constant is obtained by fitting the data 
for 10 < t < 20 ps, cf. Figure [5|d). Decay constants k are found to be identical for 
trajectories exiting through either channel (0.093 ps^^ for channel 1, Fig ^p, 0.091 
ps^^ for channel 2, FigjTjo). Again, equality of decay rates is expected on symmetry 
grounds for a symmetric ensemble. More generally, the equality of decay rates for 
trajectory subensembles reacting via distinct channels is an implicit assumption of 
statistical theories^^. An informative discussion of this equality is given in Sec. V(d) 
of ref . M 

For a system such as HON with 2 identical transition states the total forward decay 
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constant k = kf. 
B. Statistical and modified statistical rates 

Having characterized the behavior of the gap time and hfetime distributions on various 
timescales, we now consider the relation of the numerically determined decay rate to various 
statistical estimates. 

Figure [s] shows the survival probability -Ps(t) for phase points in the reactant region of 
phase space (cf. Fig. 3b of ref. I117p . This quantity is computed by Monte Carlo sampling 
phase points x in the reactant region, x G S^jflA^A, and propagating trajectories until they 
either react through either channel or a cutoff time is reached. The fraction of phase points 
surviving until time t is Ps{t)- It can be seen that Ps{t) appears to converge relatively slowly 
to a constant value /t = Ps{oo) ^ 0.9l'^. This shows that over 90% of phase points are 
trapped in the reactant region, and that the density of states pp^ = {l — /t)pa associated 
with crossing trajectories is approximately 10% of the full reactant density of states pA- 

The value of p^ has also been computed using the relation 

pI = 2s<P{E), (3.4) 

and shown to be identical with the result obtained via the survival probabilitj^J^. 

Numerical values of relevant quantities are s = 0.163 ps, (p = 0.00085 x h"^, where h is 
Planck's constant, pa = 0.795 x h^/eV and p^ = 0.0715 x h^/eV. Computed values for the 
total statistical decay rate coefficient (in the symmetric case, equal to the single channel rate 
constant kf) are kj^^^ = 0.517 ps~^ (RRKM, uncorrrected) and /c^^™ = 1/s = 5.75 ps~^ 
(RRKM, corrected). These values are to be compared with the long time decay rate k = 
0.092 ps"i . 

It is immediately apparent that, as noted previously (ref. [631 although cf. ref. [TO]), the 
value of the statistical rate constant "corrected" for the volume of trapped reactant phase 
points is both much larger than the uncorrected statistical rate and in significantly greater 
disagreement with the exact (numerical) value of k. In fact, the ratio kj'^'^^/k = 5.62, 
so that the uncorrected statistical rate coefficient is within a factor of 6 of the numerical 
escape rate. This is presumably due to compensating errors in the statistical calculation!^. 
The presence of trapped regions of reactive phase space decreases the volume of phase space 
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that is available for reactive trajectories to explore; this effect tends to increase the value 
of the actual escape rate with respect to the statistical estimate. As can be seen from the 
numerical values given above, if this were the only factor affecting the rate then the actual 
rate would be a factor of 10 larger than the RRKM estimate. However, additional dynamical 
trapping mechanisms, that are as yet not fully understood for multimode system^^^'^^'^, 
serve to delay the exit of phase points from the reactant region. The competition between 
these two effects then results in a value of the numerical escape rate that is fairly close to 
the simple RRKM estimate. 

C. Statistical branching ratio accompanied by nonexponential decay 

In Figure [9] we plot log[A^(t)] versus t for each subsensemble, where N{t) is the number 
of phase points remaining in the well at time t. Except at extremely short times ^0.5 ps, 
the two curves are essentially identical. For those times during which the decay curves are 
identical, the associated product branching ratio is equal to its statistical value, unity. The 
fact that the curves can be overlaid even at times for which the decay of each ensemble 
is nonexponential (algebraic) then implies that we have a statistical product ratio in the 
absence of exponential decay. (Note that exponential decay is usually taken to be the 
signature of "statistical" dynamic The identity of the decay curves also implies the 
equality of the algebraic decay exponents noted previously. 

A qualitative explanation of this behavior is as follows: note that we start with a highly 
asymmetric initial ensemble, entering via DSi^in only. This bundle of trajectories passes 
through the HCN well, over to DS2. Some fraction of the ensemble exits through DS2,out 
and is lost. The rest of the trajectories then turn back, and pass through the well again, 
over to DSi. Some fraction is again lost. (These are the "pulses" seen in the short-time gap 
time distribution.) And so on. 

The point is that, after the trajectory bundle has oscillated back and forth in the HCN well 
several times, the set of phase points still in the well behaves as if as if it were an ensemble 
consisting of trajectories initiated in equal numbers on both entry dividing surfaces and 
then propagated for t > 0.5 ps. The underlying idea here is that, when trajectories are 
"turned back" from a DS, their subsequent time evolution can be qualitatively similar to 
that of trajectories actually initiated on the same DS. For example, trajectories turned back 
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at DS2 that are just "outside" the incoming reactive cyhnder manifold associated with the 
dividing surfac ^ * ^^° * ^^° l will track (shado'wff^ trajectories that actually enter the reactant 
region through DS2,in inside the reactive cylinder and lie close to the boundary (unstable 
manifold W^"). 

The HCN gap time/lifetime distributions therefore imply that the initial ensemble effec- 
tively attains the full symmetry of the reactant phase space some time before it fully relaxes 
to the stage where exponential decay is observecP^, so that the system exhibits statistical 
branching ratios in the absence of exponential decay. 

D. Comparison with other calculations 

The trajectory calculations analyzed here were carried out using the MCH potential 
surfac#^ at a constant energy E = —11.88 eV (0.2 eV above the saddle energy). Tang, Jang, 
Zhao and Rice (TiZR)^! have carried out classical trajectory calculations using the same 
potential surface for a number of energies, and have applied an approximate versiorp21 of the 
3-state statistical theory of Gray and Ricd^as well as a reaction path approach®^ to compute 
isomerization rate constants for the 3 DoF, zero angular momentum HCN isomerization. 

The energy value used by TJZR closest to that of the present work is E = —11.5 eV. At 
this energy, TJZR extract an isomerization rate of 0.146 x 10~^ a.u. from their trajectory 
calculations, corresponding to a mean lifetime of 0.166 ps^^^. This rate, which determines 
the timescale for decay of the CNH population and is therefore strictly the relaxation rate 
k = kf + fc J60 | 6i | 70 |^ jg actually obtained by computing the average escape time (lifetime) 
for an ensemble of trajectories in the CNH well, while somewhat arbitrarily omitting from 
consideration short-lived trajectories with lifetimes < 1500 a.u.l^^. 

The value of the mean lifetime computed by TJZR at -E = —11.5 is close to the mean gap 
time (0.174 ps) computed at E = —11.88 eV^. Both approximate rate theories applied 
to the problem by TJZR give isomerization rates within factors of 2-3 of the trajectory 
values. These approximate statistical theories are however only capable of describing short 
time kineticJ^^'^ the slow decays apparent from the trajectory data of TJZR (see Figure 7 
oP^ cf. Figs [sj |6] and [t] of the present paper) are not predicted by 2- or 3-state statistical 
models. The algebraic decay observed at intermediate times in the present work suggests 
that incorporation of a small number of approximate intramolecular bottlenecks into the 
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statistical uiodeW^^^^^^ is unlikely to lead to an accurate description of decay dynamics 
at longer times. 

E. Relation to nonstatistical behavior in reactions of organic molecules 

Finally, it is very interesting to note that the dynamical behavior associated with the 
gap time distribution analyzed here for HCN isomerization provides an exemplary instance 
of the highly nonstatistical dynamics discussed by Carpenter in the context of reactions of 
organic molecule^^^'^. 

Thus, as discussed above, an ensemble of trajectories is launched into a reactant region 
(the HCN well) for which there are two equivalent exits ("products"). Any statistical theory 
predicts a 1:1 branching ratio, as there is by symmetry equal probability of leaving via exit 

1 or 2. The trajectory calculations however show that a significant fraction of the initial 
ensemble of trajectories simply passes through the HCN well and exits directly via channel 

2 in a very short time; the remainder of the ensemble then sloshes back and forth in the 
well leading to the pulses seen at short times in the gap time distribution; these pulses 
are associated with bursts of exiting trajectories alternating between channels 2 and 1. At 
longer times, t > 0.5 ps, we enter a regime in which decay into either channel is equally 
likely. By this time, however, most of the trajectories have already exited the well. 
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IV. SUMMARY AND CONCLUSION 



In this paper we have revisited the general approach to classical unimolecular reaction 
rates due to Thiel^^ in light of recent advances in the phase space formulation of transition 
state theory for multidimensional system^™- ^ ' ^^" ' ^^^ ' "''^ ' ^^" ' ^^^ ' "^ ' "'^ ' ^^" ! We showed that 
key concepts in Thieles's approach, namely the phase space dividing surface separating 
reactants from products, the average gap time, and the volume of phase space associated 
with reactive (interesting crossin^^ trajectories, are both rigorously defined and readily 
computed within the phase space approadJ^SHninillSQl 

The distribution of gap times is a central element of Thiele's approachP^'^. Here, we 
have analyzed in detail the gap time distribution and associated reactant lifetime distribu- 
tion for the isomerization reaction HCN ^ CNH, previously studied using the methods of 
phase space transition state theory in refs I116I117I Both algebraic (power law) and expo- 
nential decay regimes have been identified. The dynamical origins of power law behavior or 
'fractional dynamics' in multimode Hamiltonian systems, especially at intermediate times, 
remain obscur#22|_ When combined with the normal form algorithm for computation and 
sampling of the phase space dividing surface, the cumulative gap time distribution is never- 
theless a powerful diagnostic for reactive dynamics in multi-dimensional systems. 

We have also compared statistical estimates of the isomerization rat d^^ * ^^ * ^^ * ^^^ * ^^^ *with the 
numerically determined decay rate. We have found that, as noted by other^, correcting the 
RRKM estimate to account for the measure of the reactant phase space region occupied by 
trapped trajectories results in a drastic overestimate of the isomerization rate. Compensating 
but as yet poorly understood trapping mechanisms in the reactant region serve to slow the 
escape rate sufficiently that the uncorrected RRKM estimate turns out to be reasonably 
accurate, at least at the particular energy studied. 

In the planar model of HCN isomerization studied here, trajectories can exit the HCN 
well through either of 2 channels, where the channels are related by a reflection symmetry. 
Analysis of the decay properties of the subsensembles of trajectories that exit through par- 
ticular channels shows that, despite a highly asymmetric distribution of initial conditions, 
the complete trajectory ensemble effectively attains the full symmetry of the system phase 
space on a short timescale t < 0.5 ps, after which the product branching ratio is 1:1, the "sta- 
tistical" value. However, at intermediate times, this statistical product ratio is accompanied 
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by nonexponential (algebraic) decay. 

We have also pointed out the close parallels between the dynamical behavior inferred 
from the gap time distribution for HCN and nonstatistical dynamics in reactions of organic 
molecules discussed by Carpentei* ^^ * ^^ l 
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APPENDIX A: DIVISION OF REACTANT PHASE SPACE INTO TRAPPED 
AND REACTIVE COMPONENTS 

The Poincare recurrence theoremP^HSl^ used in conjunction with the geometrical struc- 
tures constructed in the phase space reaction rate theory described in Section [Tl| allows us to 
give a precise treatment of the division of the reactant phase space into reactive ( "crossing" ) 
and trapped regions. 

Suppose that, for a fixed energy E, there are d saddles associated with dividing surfaces 
DSj(£'), i = 1,. . . ,d, where each dividing surface is divided into two 2n — 2 dimensional 
hemispheres DSj^in(£') and DSj^out(-E') that control entrance and exit to a compact reactant 
region. We first show that every trajectory (except for a set of measure zero) that enters 
the reactant region through a dividing surface DSj^in(-E') will exit the reactant region at a 
later time. 

To this end, consider a set V of {q,p) in DSi,in(-E') of positive volume with respect to 
the (Lebesgue) measure dcr = d"~^gd"^^p. Suppose that the points in y as a set of initial 
conditions for Hamilton's equation give trajectories which stay in the reactants region for all 
time ip > 0. Then the region swept out by these trajectories has infinite volume with respect 



to the (Lebesgue) measure da A dip (see Sec. II B). This contradicts the compactness of the 



reactant region. We can thus conclude that every initial condition on DSj.in(-E) (except 
for a set of measure zero with respect to the measure da = d"~^gd"'~^p) gives a trajectory 
which leaves the reactant region at a later time. This conclusion holds for both compact 
and noncompact (dissociative) product regions. 
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In the case that the product regions are also compact, we can invoke the Poincare recur- 
rence theorem for Hamiltonian dynamics on compact energy surface#^ Consider any open 
set in a compact energy surface. Then, with the possible exception of a set of (Lehesgue) 
measure zero, trajectories of Hamilton's equations with initial conditions starting in this open 
set return infinitely often to this set. 

Consider DSj4n(-£'), for any i. Trajectories starting on this surface must enter the reactant 
region. That is, the vector field defined by Hamilton's equations, evaluated on DSj in(-E), 
is transverse to DSj^in(-E) and pointing strictly into the reactant region (as proved irjn^. 
This is the mathematical manifestation of the "no-recrossing" property. Since the vector 
field defined by Hamilton's equations points strictly into the reactant region on DSj^nl-^') 
we can construct a "thin" open set, Oi-\-a{E)^ containing DSi,in(-E') having the property that 
all trajectories starting in this open set also enter the region. 

Now, by construction, no trajectory leaving (9j.in(-E') and entering the reactant region can 
ever intersect (9j.in(i?) without first leaving the region (which must occur through DSj^out (-£'), 
for some j). The reason for this is that the vector field defined by Hamilton's equations 
restricted to (9j;in(-E) is pointing strictly into the reactant region. By the Poincare recurrence 
theorem, with the possible exception of a set of zero Lebesgue measure, every trajectory 
starting on 0^\a{E) intersects Oi-^i^{E) infinitely often. Therefore, we conclude that, with 
the possible exception of a set of zero Lebesgue measure, every trajectory that enters the 
reactant region exits and re-enters the region an infinite number of times. We can summarize 
as follows: Almost all trajectories that enter a given reactant region of phase space exit the 
region at a later time. Moreover, after exiting, they will re-enter the same region at a later 
time, and this "entrance- exit" behaviour continues for all time thereafter. 

We can state this result also in a slightly different, but equivalent, way: Almost all 
trajectories that exit given reactant region will return to the same region at a later time. 
Moreover, after returning, they will exit the region again at a later time, and this "exit- 
return" behaviour continues for all time thereafter. 

The immediate implication is that, with the possible exception of a set of Lebesgue 
measure zero, no trajectory can escape the reactant region of phase space that is not in the 
reactive (crossing) volume. The further implication is that, with the possible exception of 
a set of Lebesgue measure zero, the volume of the reactant region of phase space consists 
of two components- the reactive volume and the trapped volume. The boundary between 
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these two sets is hable to be exceedingly complicated (fractaP^. 

APPENDIX B: DERIVATION OF THE DELAY DIFFERENTIAL EQUATION 
FOR 



In this Appendix we give a concise derivation of the delay differential equation satisfied by 
the reactant lifetime distribution P(t) in the DLGM of Dumont and Brumei'^. We consider 
a fixed value of the energy E. 



In the notation of Sec. |II C[ the fraction of trajectories initiated on DSin with gap time 
s>t' is 

r+oo 

/ dsV{s) = sF{t'), (Bl) 
Jt' 

so that the probability of a trajectory initiated on DSin having a gap time t + t' < s < 
t + t' + dt, t > 0, given that s > t' is 

s P(t') ^ ' 

The condition for statistical decay given by Dumont and Brumer can then be written!^ 

^^^^ = P(t), Vt > 0, > r (B3) 

where r is the relaxation time. The meaning of this condition is as follows: consider those 
trajectories on DSm with gap time s > t'; the fraction of these trajectories having gap times 
t + t' < s < t + t' + dt is equal to P(t)dt, the fraction of reactant phase points with lifetimes 
t t + dt, for t' > T. An equivalent formulation of the condition for statistical decay is 

F{t + t') = -P(t)P(t') \/t >0,t' > T. (B4) 



Conditions (B3) and (B4) are clearly satisfied in the case of an exponential lifetime distri- 
bution, P(t) = ke""'*. 

To obtain the Delayed Lifetime Gap ModeP^ for the lifetime distribution: 

• Consider the gap time distribution V{t) for the statistical component of reactive phase 
space only. The reactive phase space has to be partitioned into a direct and a statistical 
component, perhaps using a trajectory divergence criteriorpSl. 
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• Set ks = (s) ^, where the mean gap time is evaluated by averaging over the statistical 
component only. We therefore have P(0) = ks 

• Assume that there are no gaps in the statistical component for s < t. That is, 'P(s) = 
for s < r which implies F{t) = k^, < t < t. 

• Set t' = T (the relaxation time) to obtain the delay differential equation for the DLGM 
lifetime distribution: 

^ P(t + r) = -A;s P(t) Vt > (B5) 

CJ.6 

This equation can be solved using the Laplace-Fourier transformP^. 
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FIGURE CAPTIONS 



FIG. 1: Phase space structures for unimolecular reaction (schematic), (a) Definition of 
reactant region, NHIM and dividing surface DS(i?) = DSin(-E') UDSout(-£')- (b) Definition of 
gap time s and hfetime t. 

FIG. 2: Isopotential surfaces of the HCN potential energy surface of ref. llSUl in polar repre- 
sentation of the Jacobi coordinates r, R, and 7. 

FIG. 3: HCN gap time distribution V{s). (a) Complete ensemble, (b) Subensemble reacting 
via channel 1. (c) Subensemble reacting via channel 2. 

FIG. 4: HCN cumulative gap time (reactant lifetime) distribution F{t) at short times. 
LogF(t) is plotted vs t for < t < 0.4 ps. (a) Total ensemble, (b) Subensemble reacting 
via channel 1. (c) Subensemble reacting via channel 2. 

FIG. 5: Cumulative gap time (reactant lifetime) distribution F{t) for the complete ensem- 
ble, (a) A log-log plot shows power law decay at intermediate times, (b) Log plot shows 
exponential decay (10 < t < 20 ps). 

FIG. 6: Cumulative gap time (reactant lifetime) distribution F{t) for the subensemble 
reacting via channel 1. (a) A log-log plot shows power law decay at intermediate times, (b) 
Log plot shows exponential decay (10 < t < 20 ps). 

FIG. 7: Cumulative gap time (reactant lifetime) distribution F{t) for the subensemble 
reacting via channel 2. (a) A log-log plot shows power law decay at intermediate times, (b) 
Log plot shows exponential decay (10 < t < 20 ps). 

FIG. 8: HCN survival probability -Ps(t). Psit) is the fraction of an ensemble of trajectories 
uniformly distributed throughout the HCN region of phase space at t = remaining in the 
well at time t. Trajectories are removed from the ensemble once they exit the HCN region 
by crossing DSj^out, j = 1, 2; they cannot re-enter the region. 

FIG. 9: The log of the number N{t) of trajectories remaining at time t versus t is plotted 
for each subensemble: channel 1 (blue), channel 2 (red). 
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